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We investigate the bifurcation cascades of a linear librational orbit in a generalized 
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5— i . class of Henon-Heiles potentials. The stability traces of the new orbits born at its 
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<^ I bifurcations are found numerically to intersect linearly at the saddle energy (e = 1), 

forming what we term the "Henon-Heiles fans". In the limit close to the saddle 
Q ■ energy (e — > 1), where the dynamics is nearly chaotic, we derive analytical asymptotic 
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expressions for the stability traces of both types of orbits and confirm the numerically 
determined properties of the generalized "Henon-Heiles fans". As a bonus of our 
results, we obtain analytical approximations for the bifurcation energies e n which 



become asymptotically exact for e n — > 1. 
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I. INTRODUCTION 

The approximation of the exact density of states of a quantum system in terms of classical 
periodic orbits via semiclassical trace formulae is a fascinating subject which has triggered a lot 
of research (see [1, 2] and the literature quoted therein). It presents a nice illustration of the 
correspondence between classical and quantum mechanics, besides allowing one to approximately 
determine quantum shell structures in terms of classical mechanics (see [2] for applications in 
various fields of physics) . In Hamiltonian systems which are classically neither regular nor purely 
chaotic, this semiclassical theory is enriched - but also complicated - by the many facets of non- 
linear dynamics. One of them is the bifurcation of periodic orbits when they undergo changes of 
stability [3]. 

An essential ingredient to determine the stability of a periodic orbit is its so-called stability 
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matrix Mj_, appearing in the amplitudes of Gutzwiller's trace formula [4], which is determined 
from the linearized equations of motion around the periodic orbit. The analytical calculation of 
Mi for non-integrable systems with mixed dynamics is in general not possible; the only non-trivial 
example is, to our knowledge, that of a two-dimensional quartic oscillator [5]. 

In this paper we investigate the stability matrix Mj_ of the simplest orbit in a class of two- 
dimensional potentials which are a generalization of the famous Henon-Heiles (HH) potential [6] 
that has become a text-book example of a system with mixed classical dynamics. For small 
energies the motion is dominated by a harmonic-oscillator part and is quasi-regular; at energies 
close to and above the saddles (e = 1), over which a particle can escape, the motion is quasi- 
chaotic (see, e.g., [1, 2, 6, 7] and the literature quoted therein). At all energies below the saddle, 
there exists a straight-line librating orbit A oscillating towards one of the saddles. This orbit 
undergoes an infinite sequence of stability oscillations and hence a cascade of bifurcations, which 
can be understood as the main mechanism of the transition from regular motion to chaos [8-10]. 
The stability traces of the new orbits R and L generated at the bifurcations are found numerically 
[9] to intersect linearly at the saddle energy (e = 1), forming what was termed the "Henon-Heiles 
fans" [11]. 

In the present paper we present analytical calculations of the stability traces of both the A 
orbit and the new orbits R and L bifurcating from it. The results are obtained in the limit close 
to the saddle (e — > 1) and hence asymptotically valid as the bifurcation energies e n approach the 
saddle energy e — 1. They confirm analytically the numerical properties of the "Henon-Heiles 
fans" also in the generalized HH potential. As a bonus, we obtain analytical expressions for the 
bifurcation energies e„, which are mathematically valid asymptotically for e n — > 1, i.e., for n — > oo, 
and practically for n > 7 within 5 digits. 

In Sec. II we present the generalized Henon-Heiles system and discuss its shortest orbits, the 
bifurcation cascade of the linear A orbit and, in particular, the properties of the "Henon-Heiles 
fans". In Sec. Ill we present the basic ideas of our analytical approach and the essential results, 
while the technical details of our calculations are given in the Appendices A and B. In Sec. IV we 
present an alternative perturbative approach for evaluating the stability traces, with the details 
given in Appendix C, and compare its results with those of the non-perturbative calculations. 
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II. BIFURCATION CASCADES IN THE HENON-HEILES SYSTEM 
A. The generalized Henon-Heiles Hamiltonian 

In this paper we investigate the following family of Hamiltonians: 

(1) 



H G hh = \ (pl+p 2 y) + \ {x 2 + y 2 ) + a 



~y 3 + ix 2 y 



where 7 > is a parameter specifying specific members of the family, and a > is a chaoticity 
parameter that can be scaled away with the energy as shown below. For 7=1, the Hamiltonian 
(1) reduces to the standard Henon-Heiles (HH) Hamiltonian [6]; we therefore call (1) here the 
"generalized Henon-Heiles" (GHH) Hamiltonian. The HH system with 7 = 1 has symmetry: 
it is invariant under rotations around the origin by 27r/3 and 47r/3, and under reflections at three 
symmetry lines with the angles ±7r/6 and n/2 with respect to the x axis. It exhibits three saddles 
at energy E sad = l/6a 2 , the equipotential lines at E = E sad forming an equilateral triangle. For 
7^1, the symmetry is lost and only the reflection symmetry at the y axis remains; there 
are, however, still three saddles over which the particle can escape. For 7 = the system becomes 
separable and has only one saddle on the y axis (cf. [12, 13]). 

After multiplying the Hamiltonian (1) by a factor 6a 2 and introducing the scaled variables 
x',y\e by 

x' = ax, y' = ay, e = 6a 2 E = E/E sad , (2) 

the scaled Hamiltonian becomes independent of a, and for a given 7 there is only one parameter 
e that regulates the classical dynamics. For simplicity of notation, we omit in the following the 
primes of the scaled coordinates x, y but keep using the scaled energy e. 

For 7 = 1, the three saddles are at the scaled energy e = 1; one of them is positioned at x — 0, 
y — 1. For 7^1, the saddle with energy e = 1 persists at the same position, while the two 
other saddles lie at different energies and are positioned symmetrically to the y axis. For a more 
detailed description of the topology of the potential (1) (and an even larger class of generalized 
HH potentials) and its shortest periodic orbits, we refer to a forthcoming publication [14]. 

The shortest periodic orbits of the standard HH system (7 = 1) have been extensively discussed 
in the literature [8-10, 15], and their use in semiclassical trace formulae for the quantum density 
of state of the HH system was investigated in [12, 16-18]. 
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B. The motion along the A orbit 

As mentioned above, we use henceforth the symbols x, y for the scaled coordinates (corre- 
sponding to a — 1), along with the scaled energy e given in (2). The equations of motion for the 
Hamiltonian (1) are then 

x(t) + [l + 2>yy(t)]x(t) = 0, (3) 
y(t)+y(t)-y 2 (t)+ 1 x 2 (t) = 0. (4) 

In the present work we focus on the linear orbit that librates along the y axis, here called the A 
orbit. It goes through the origin (x,y) = (0,0) and towards the saddle at (x,y) = (0,1) which 
it, however, only reaches asymptotically for e — > 1 with a period Ta — > oo. Since this orbit has 
x A (t) = XA(t) = at all times t, its equation of motion is 

y A (t)+y A (t)-y A (t) = 0, (5) 

which can be solved analytically [12]. We give here the result in the most general form, relevant 
for our following development, where the initial point along the y axis is given as y = y A (t = 0). 
The solution is then: 

y A (t) =yi + (y2-yi)sn 2 (z,K) , (6) 

z = a K t + F((p, k) . (7) 

Here sn(z, k) is a Jacobi elliptic function [19] with argument z; its modulus k and the constant a K 

are given by 

k=J— — — , a K = ^/{y3-yi)/Q, (8) 
V I/3 — I/1 

in terms of the three real solutions of the equation e = 3y 2 — 2y 3 < 1 given by 

Vl = 1/2 - cos(tt/3 - 0/3) , y 2 = 1/2 - cos(tt/3 + 0/3) , y 3 = 1/2 + cos(0/3) , (9) 

with cos0 = 1 — 2e. The function F(<p,k) in (6) is the incomplete elliptic integral of first kind 
with modulus k, the argument ip being determined by the initial condition: 



. yo-yi nn x 

09 = arcsin J . (10) 

V 2/2 - y\ 

yi and y 2 are the lower and upper turning points, respectively, of the A orbit along the y axis. 
The period and the action of the (primitive) A orbit are given by 

T A (e) = — K(«) , S A (e) = ^ [E(«) + c K K(«)] , (11) 
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with c K = — 2(y 3 — y 2 )(2y3 — y 2 — 2/i)/9, in terms of the complete elliptic integrals of first and 
second kind, K(k) and E(/t) (we use the notation of [19]). 

Note that in the limit e — * 1, we have y 2 — ► 1, 2/3 — ► 1 and « — > 1, so that K(/c) and diverge 
(while Sa remains finite). The A orbit then is no longer periodic (and may be called a "homoclinic 
orbit" [3]). Expanding Ta around e = 1, one finds the asymptotic form [9] 



C. The bifurcation cascade of the A orbit in the standard HH potential 

While approaching the saddle as e — > 1, the A orbit undergoes an infinite cascade of pitchfork 
bifurcations, giving birth to a sequence of new orbits R 5 , L 6 , R7, L 8 , .... This scenario, which 
has some similarities to the Feigenbaum scenario [20], was discussed extensively in [9], and the 
analytical forms of the newborn R and L orbits in terms of periodic Lame functions were discussed 



In Fig. 1 we show the traces of the stability matrix Mj_, defined in (14) below, of the A orbit 
and the orbits bifurcated from it, plotted versus energy e. Whenever trMj_ = 2, a bifurcation 
occurs. We see the successive bifurcations at increasing energies e n ; upon repeated zooming the 
upper end of the energy scale near e = 1 (from bottom to top), the pattern repeats itself in a self- 
similar manner. The bifurcation energies e„ form a geometrically progressing series (see [9, 10] for 
details), cumulating at the saddle energy (e = 1) such that es(Rs) < ee(L 6 ) < e^iRr) < • • • < 1, 
where the parentheses contain the names of the new orbits born at the pitchfork bifurcations. 
These are alternatively of R type (rotations) and of L type (librations). (The subscripts in the 
orbit names indicate the Maslov indices appearing in the semiclassical trace formulae; the index 
of the A orbit increases by one unit at each bifurcation.) Due to the discrete symmetries of the 
system, all these pitchfork bifurcations are isochronous and hence not generic (cf. [14]). 

In Fig. 2, we show again trMj_ - in the following briefly termed the "stability traces" - of the 
same orbits, but this time plotted versus their respective periods T. On this scale, tr Mj^Ta) 
(shown by the heavy line) is numerically found [8] for large Ta to go like a sine function; its period 
AT = 3.6276 was shown in [9] to be given analytically by AT = 2n A/3. The exact calculation 
of the function tr Mj_a(Ta) is, however, not trivial at all. It is one of the objects of our present 
investigations (see Sec. Ill B) . 





1) 



(12) 



in [10]. 
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D. The "Henon-Heiles fans" 

An interesting property of the stability traces of the R and L orbits born at the bifurcations, 
which has been observed numerically [9] and termed the "Henon-Heiles fan" structure [11], is 
emphasized in Fig. 3. Here we plot the stability traces of the primitive A orbit and the first three 
primitive pairs of R and L orbits versus the scaled energy e. We note two prominent features 
(which can also be recognized in Fig. 1): 

(i) The functions trM^ RL (e) are approximately linear up to (and even beyond) the barrier energy 
e = I. 

(ii) The curves tr M ifl) i(e) intersect at e = 1 in one point each for all R and L type orbits with 
Maslov indices > 8, positioned at the values 2 ± d with d = 6.183 ± 0.001, thus forming two 
fans emanating from these points. The uncertainty in the parameter d comes from the numerical 
difficulty of finding periodic orbits (which was done using a Newton-Raphson iteration procedure) 
close to bifurcations; our result for d was obtained for R„ and L„/ orbits with 9 < n, n! < 13, 
evaluated at e = 1. The upper limit n = 13 is due to the numerical problems only; we expect that 
the same value d = 6.183 ± 0.001 holds also for all higher n. 

We found exactly the same types of "HH fans" for the generalized HH systems given by the 
Hamiltonian (1) for the bifurcation cascade of the A orbit along the y axis, whereby the slopes of 
the fans and hence the value of d depend on the parameter 7. The "GHH fans" can be described, 
for large enough n, by the empirical formula 

trMSg(e) = 2 =F c RL (l) j^j , (e > e n ) (13) 

where the negative and positive sign belongs to the R and L type orbits, respectively. At e = 1 
the curves tr MxR^e) intersect linearly at the two values trMjj^l) = 2 =1= crl(i), so that the 
parameter d given above for the standard HH potential is d = crl(1). The numerical values for 
c rl(i) are shown by crosses in Fig. 5 below. 

The main goal of our paper is to find analytical support for these numerical findings. In Sec. 
Ill we will, indeed, confirm the empirical formula (13) analytically in the asymptotic limit e — > l. 

III. ASYMPTOTIC EVALUATION OF STABILITY TRACES 

In this section we derive analytic expressions for trMj_(e) of the A, R and L orbits in the GHH 
system, which are valid in the asymptotic limit e — > 1, i.e., close to the barrier. Before presenting 
them in Sec. Ill B and Sec. Ill C, we recall the definitions of the stability matrix trMj_ and of the 
monodromy matrix M of which it is a submatrix. 



A. Monodromy and stability matrices 
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1. Stability matrix and the Hill equation 

The analytical calculation of the stability matrix Mi of a periodic orbit in a non-integrable sys- 
tem is in general a difficult task. We recall that the stability matrix is obtained from a linearization 
of the equations of motion and defined by 

<^(T) = M ± 5£ ± (0), (14) 

where 6£±(t) is the (2N — 2)-dimensional phase-space vector of infinitesimally small variations 
transverse to the given periodic orbit (N being the number of independent degrees of freedom), 
and T is the period of the orbit. For N = 2 dimensional systems, we may choose £±(t) = {q,p) 
where q is the coordinate and p the canonical momentum transverse to the orbit in the plane of 
its motion. (q,p) then form a "natural" canonical pair of Poincare variables, normalized such that 
(q,p) = (0,0) is the fixed point of the periodic orbit on the projected Poincare surface of section 
(PSS) . For two-dimensional Hamiltonians of the form "kinetic + potential energy": H = T + V 
(and particles with mass m — 1, so that p — q), the Newtonian form of the linearized equation of 
motion for q(t) becomes the Hill equation (see the text book [21] for an explicit discussion) 

q(t) + V qq (t)q(t) = 0, (15) 

where V qq (t) is the second partial derivative of the potential V with respect to q, taken along the 
periodic orbit, and the two-dimensional stability matrix is given by 

(16) 

For isolated periodic orbits, solutions of (15) with q(t) ^ are in general not periodic. However, 
when the orbit undergoes a bifurcation, (15) has at least one periodic solution which describes 
the transverse motion of the new orbit born at the bifurcation; the criterion for the bifurcation to 
occur is trM_L = +2 (cf. [21]). 

For particular systems, the Hill equation (15) may become a differential equation with known 
periodic solutions. For the GHH systems under investigation here, the Hill equation for the A 
orbit directed along the y axis is given by (3), with y(t) replaced by in (6), and becomes the 

Lame equation (see, e.g., [22]) whose periodic solutions are the periodic Lame functions (see [10] 
for the details). However, the elements of Mj_ in (16) can in general not be found analytically. One 
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of the rare exceptions is that of the coupled two-dimensional quartic oscillator for which Yoshida 
[5] derived an analytical expression for trMj_ as a function of the chaoticity parameter (cf. [23]). 

Magnus and Winkler [21] have given an iteration scheme for the computation of trMj_ for 
periodic orbits in smooth Hamiltonians. We have tried their method for the A orbit in the HH 
system, but we found [24] that its convergence is too slow for computing tr M±A(e) with a sufficient 
accuracy that would allow to deduce the properties of the HH fans. However, in the limit e — > 1, 
it is possible to use an asymptotic expansion of the function sn appearing in yA(t) of (6), which 
allows us to compute trM^(e) analytically, as discussed in Sec. Ill B. 



2. Matrizant and monodromy matrix 

For curved periodic orbits - such as the R and L type orbits bifurcating from the A orbit 
in the GHH systems - which usually can only be found numerically, the phase-space variables 

transverse to the orbit used in the definition (14) of the stability matrix Mi cannot be con- 
structed analytically. Instead, one must in general use Cartesian coordinates and resort to the full 
monodromy matrix M defined below. For N = 2, one first linearizes the equations of motion to 
find the matrizant X(t) which propagates small perturbations of the full phase-space vector £(£) 
defined by 

£(0 = {x(t),y(t),x(t),y(t)} (17) 
from their initial values at t = to a finite time t: 

5£(t) = X(t) 5m . (18) 

For a Hamiltonian of the form H(x,y,x,y) = \ (x 1 + y 2 ) + V(x, y), the differential equation for 
X(t) is 

^X(f) = I - |X(/) (19) 




with the initial conditions 



X(0)=I 4 , (20) 

where I2, I4 are the two- and four-dimensional unit matrices and U(t) is the two-dimensional 

Hessian matrix of the potential, taken along the periodic orbit (po): 

d 2 V 

Uij(t) = Q X .Q X . { x ( t )>y( t )}p° ' ( x *> x 3 = x ^y)- ( 21 ) 

Having solved (19), the monodromy matrix M of the given periodic orbit with period T is 
defined by 

M = X(T) . (22) 
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In an autonomous system, M has always two unit eigenvalues corresponding small initial variations 
along the periodic orbit and transverse to the energy shell. After a transformation to an "intrinsic" 
coordinate system, in which one of the coordinates is always in the direction r\\ (with momentum 
p\\ = f\\) of the periodic orbit [4], M can be brought into the form 

/ Ml ... \ 

' 1 ... 
1 



M 







V 



(23) 



/ 



where the dots denote arbitrary non-zero real numbers and Mi is the stability matrix. The 
diagonal elements in the lower right block of (23) then correspond to 



(24) 



The transformation to such an intrinsic coordinate system is quite non-trivial [25] and not 
unique. For curved orbits it can in general only be found numerically and is therefore not suitable 
for analytical calculations. For the curved R and L orbits of our system, we therefore have to 
resort to the full monodromy matrix M (22) via the solution of (19). For the evaluation of their 
stability traces, we only need the diagonal elements of M and can then use the obvious relation 
trM, =trM-2. 



B. Asymptotic evaluation of the stability trace trMj^(e) for e — > 1 

In the limit e — > 1, where the modulus re defined in (8) goes to unity, we may approximate 
DA(t) by the leading term in the expansion of the function sn(z, re) around re — 1 (see [19]) 

sn(z, re) k, tanh(z) . (re — > 1) (25) 

Since the function tanh(z) is not periodic, we have to approximate yA(t) in two portions. Taking 
t 2 as the time where the orbit passes through its maximum at y 2 , i-e., 

VA(t 2 ) =V2 h = [K(re) - F(<p, K)]/a K , (26) 

we define the asymptotic expression for the A orbit over one period by 

y A (t) = G(t 2 - t) Y^t) + G(t - t 2 ) Y 2 (t) , < t < T A , (27) 

where the functions Yi(t) and Y 2 (t) are given by 



Y^t) = Vl + (y 2 - yi) tanh 2 (2;) , Y 2 {t) = Vl + (y 2 - yi) t^h\z - 2 K(re)) , (28) 
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with z given in (7). Although the function (27) is not analytic at t — t2, it suffices to find an 
asymptotic expression for trM± A (e) valid for e — > 1. 

The details of our calculation are given in Appendix A. The analytical asymptotic result for 
trMxA(e) is given in (A26) in terms of associated Legendre functions. In the limit e — > 1, the 
energy dependence of trM^ A (e) goes only through the period T\(e): 

trM ±A (e) » trMH(e)=trMH(T A (e), 7 ), (e -> 1) (29) 

where trM^ (T A , 7 ) is a universal function given by 

trM^(T A , 7 ) = +2 |F A ( 7 )| cos[v^T2^T A - $ A ( 7 )] . (30) 

The phase function $a( 7 ) is defined through Eqs. (A28) and (A30), and the amplitude function 
|F A ( 7 )| is given explicitly in (A31). We recall that 7 is the potential parameter of the GHH 
potential (1) is 7 = 1 for the standard HH potential. For this case, the result (30) becomes 

trMfJ (T A , 1) = 2.68043976 cos(VST A + 1.56782696) , (31) 

where the numerical constants have been calculated for 7 = 1. The period of the cos function 
in (31) was correctly shown in [9] to be 27rA/3, but the phase &a(i = 1) and the amplitude 
2|F A ( 7 = 1)| were only obtained numerically. The asymptotic relation (29) had already been 
observed numerically in [8, 9]. 

The result (31) is shown in Fig. 4 by the dotted line and compared to the exact numerical result 
from [9], shown by the solid line. We see that the agreement becomes nearly perfect for Ta > 10.5, 
corresponding to e > e 6 . The asymptotic result (29), (30) allows us to give analytical expressions 
for the bifurcation energies e„ in the asymptotic limit e n — > 1. The pitchfork bifurcations of the 
A orbit occur when trM^ A = +2. We therefore define approximate bifurcation energies e* by 

trMH(T A (e;), 7 ) = +2. (32) 

Using the asymptotic form of T A (e) in (12) and (30), we can give the solutions of (32) in the 
following formulae 

~ l-432 exp{-[$ A ( 7 )-arccos(l/|F A ( 7 )|) + 27rA;]/ v ^T2^}, (R) 



e* 2k ^ 1 -432 exp{-[$4( 7 ) + arccos(l/|F A ( 7 )|) + 2%k]/y/l + 2 7 }, (L) (33) 

where k — 3, 4, 5, . . . , and the odd numbers n = 2k — 1 refer to the R type and even n = 2k to the 
L type bifurcations, respectively. For e* sufficiently close to 1, i.e., for large enough n the above 
values should reproduce the numerically obtained "exact" values e n . 
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This is demonstrated for 7 = 1 in Tab. I. In the second column we give the resulting values 
of e* with 5 < n < 16 for the standard HH system, and in the third column we reproduce their 
numerical values e„ obtained in [10] as roots of the equation trM^A(e n ) = +2. As we see, the 
asymptotic results e* approach the numerical values e n very well already starting from n = 7, as 
could be expected from Fig. 4. In view of the numerical difficulties in determining the e n from a 
search of periodic orbits (cf . the remarks after Fig. 3) , the agreement is very satisfactory for all 
n > 7. 

This is in itself a remarkable result, because we are not aware of any analytical results for bifur- 
cation energies (or bifurcation values of any chaoticity parameter) in non-integrable Hamiltonian 
systems, except for the coupled two-dimensional quartic oscillator (see [10, 23]). In the present 
case, the bifurcation energies e n can be related to the eigenvalues of the Lame equation. These 
can, in principle, be given by infinite continued fractions [26], but their determination is hereby 
only possible numerically by iteration, which becomes even less accurate than the numerical so- 
lution of tr M_i_A(e n ) = +2 as done in [10]. The analytical expressions (33) therefore represent an 
important achievement of this paper. 

C. Asymptotic evaluation of trMxR i(e) for e — > 1 

For the stability traces of the R and L orbits we need, as mentioned in Sec. Ill A 2 above, to 
know the diagonal elements of the full monodromy matrix M, i.e., the elements Xa(t = T) with 
i — x,y, x, y. Since the equations (19) couple all 16 elements of X(t), this is still a considerable task. 
It can, however, be simplified considerably in the asymptotic limit e — > 1. First, we can make use 
of the "frozen y motion approximation" (in short: "frozen approximation", FA) introduced in Refs. 
[9, 10]. It exploits the fact that near the bifurcation energies e n at which the R and L orbits are 
born, their motion in the y direction is close to that of the bifurcating A orbit and, for increasing 
energy e, changes only very little. It can be shown (cf. [23] and Sec. IV below) that this may 
correspond to the first order in a perturbative expansion in the parameter e — e n , valid to leading 
order in the small quantity 1 — e n . Second, we can exploit some symmetry relations between 
the elements of M if the initial point at t — for the calculation of X(t) is chosen as the upper 
turning point in the direction of the A orbit, i.e., its maximum along the y axis. These symmetry 
relations are derived in Sec. Bl; their main consequence is that we only need to calculate the 
4x4 submatrix of Xjj(t) with spatial indices i,j = x,y, and that we have the asymptotic equality 
trM±R,L ~ 2M ra for e — > 1, see (B23). As shown below, these symmetry relations can be used 
also beyond the FA, and only in order to simplify them some properties of the FA will be exploited 
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in our further derivations. 

With these approximations, the calculation of trMj^^e) proceeds similarly as that of 
trMj_yi(e) discussed in the previous section; its details are presented in Appendix B. The an- 
alytical result is given in (B46) in terms of associated Legendre functions. After their expansion 
in the asymptotic limit e-> 1 we obtain the result 

trM^(e) = 2 =F c RL ( 7 ) j^j , (e > e n - 1) (34) 

where the "— " and "+" sign belongs to the R and L orbits, respectively. The slope function crl(i) 
is found analytically to be 

sinh[27iyi + 27J 

Eq. (34) has exactly the functional structure of the empirical "GHH fan" formula (13). Math- 
ematically, it holds asymptotically in the limit e n — > 1 to leading order in the small parameter 
VI — e n . We emphasize that this result confirms also the numerical finding that, for large enough 
n (practically, for n > 8) the functions trM^^ i(e) are linear in e from e„ up to at least e = 1. 

In Fig. 5 we show by crosses the values of crl(i), evaluated numerically from the stability of 
the R and L orbits at e = 1, as a function of 7. The solid line shows the analytical result (35). 
In the lower part of the figure, we show the region of small 7. The curve 0^(7) goes through 
zero with a finite slope which can easily be found by Taylor expanding (35) after the replacement 
cosh (71-^487 — 1/2) — > cos (flVl — 487/2). The slope at 7 = becomes 

48tt 



C -^^mh|f^2,| COSh (f^ 1 )- < 35 > 



x = 0.56320942 . (36) 
7=0 smh(27r) v ' 



This value is found analytically [11] from a semiclassical perturbative approach, in which the 
term r yx 2 y of the Hamiltonian (1) is treated as a perturbation. Using the perturbative trace 
formula given by Creagh [27] one can extract the stabilities of the R and L orbits which in this 
approach are created from the destruction of rational tori (see [11] for details). To first order in 
the perturbation, one obtains exactly the correct linear approximation to crx,( 7 ), with the slope 
(36), shown in the lower part of Fig. 5 by the dotted line [28]. 

The theoretical value of crl(1) = 6.18199717 agrees very well with the value d = 6.183 ± 0.001 
that was found from the numerical stabilities of the R n and L„/ orbits in the standard HH potential 
(7 = 1) for 9 < n, n' < 13, evaluated at e = 1. 

Our result (34) obeys a known "slope theorem" for pitchfork bifurcations [14, 29, 30]. It states 
that the slope of trMj_(e) of the new orbits born at the bifurcation point e n equals minus twice 
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that of the parent orbit. Specifically in the present system, it says 

^ tr M ±R , L (e n ) = -2 tr M ±A (e n ) . 



(37) 



We can easily obtain the slopes of trM^(e) at e = e n from the asymptotic result for trM^(e) 
given in (30) . By its Taylor expansion around the asymptotic bifurcation energy e* given by (33) , 
we find up to first order in e — e* 



trM 



(as) 
LA 1 



e) = 2 ± c a (t) 



+ o 



e - e 



(38) 



(i-<) ' ~ L(i-<) 3/2 . 

The alternating sign of the linear term is "+" for the R and "— " for the L type orbit bifurcations 
and thus opposite to that in (34). The slope function c A {"f) is found to be 
d 



ca(i) = 



trMH(T A , 7 ) 



dT A 



T A =T A (e* n ) 



= 2v^T2^VI^(7)l 2 -l 



2^/1 + 2^ 



cosh (1^487- 1) 



(39) 



sinh[27rVl + 27] 

where |i 7 U(o / ) I is given in (A31) and T A (e) in (12). Note that 0,4(7) does not depend on the 
bifurcation energy e* since trM^ (7^,7) is a periodic function of Ta- Comparing Eqs. (35) and 
(39), we see that crl{i) = 2^(7) so that the theorem (37) is, indeed, fulfilled with the correct 
sign. 



IV. PERTURB ATIVE EVALUATION OF trM ±J? L (e) NEAR e = 1 

Here we present an iterative perturbative approach for the calculation of the stability trace of 
the new orbits born at the bifurcation energies e n of the A orbit for the GHH Hamiltonian (1), 
taking R orbits as example. This approach can be useful for Hamiltonians for which we do not find 
symmetry properties like those given in (B17) and (B18), which allowed for a non-perturbative 
calculation of the stability traces. 

As the small perturbation parameter we introduce the available energy above the bifurcation 
point 

e = e - e n , (40) 

which is always positive. The x and y coordinates of the new orbits, labeled x po and y po , and 
the relevant elements of their monodromy matrices, all as functions of time t, can be expanded in 
powers of small perturbation parameter e: 

y po (t) = y A (t) +ey$(t) + ..., x po (t) = u po [xf}{t) +ex$(t) + ...], (41) 



X«(f) = Xl 0) (t) + eX^(t) + ..., X<,-(f) = u po XS'V) +eX£>(t) + ... (i * j) , (42) 
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where i, j = x,y. The superscripts ^ indicate in an obvious manner the power e m at which the 
corresponding terms appear at the m-th order of the expansion. The normalization constants u po 
of x po (t) are given by 



Note that they both are proportional to y/e, so that x po (t) goes to zero in the limit e — > e n . The 
solution of the equations (B21) with the initial conditions (B22) for the stability trace trMi^i(e) 
using the perturbative expansions (41) and (42) is presented in the Appendix C for the case of 
the R type orbits. The calculation for the L type orbits is completely analogous. The asymptotic 
result for trMj^e) is given in (CIO). 

We now compare the non-perturbative result (B42) and the perturbative approximation (CIO) 
for the stability traces trMi^i(e) with numerical results. Fig. 6 shows by solid lines the asymp- 
totic analytical results (B42) for the case 7=1. They form the "HH fans" with their linear energy 
dependence of tr M_|_ Rii (e) around e = 1, intersecting at the values tr Mj^^l) — 2 = =FCrl(1) 
with crl(1) ~ 6.182 for the R and L type orbits, respectively. As seen from this figure, they 
become approximately symmetric with respect to the line trMj_ = +2, starting from n = 9 in 
good agreement with the numerical results [10]. Note that the linear dependence of tr M±r,l(g) 
(B42) is obtained up to terms of relative order y/1 — e n . The perturbative result for the R type 
orbits (CIO) is shown by the dashed lines, in good agreement with the analytical result (B42) 
already for n > 9. 

For further comparison with numerical results, we define the "slope parameter" 



evaluating tr M± R ^ L at the barrier (e = 1) for a given orbit R n or L n born at the bifurcation energy 
e n . As shown in Sec. Ill C and Appendix B 2, this parameter tends to the asymptotic limit 0^(7), 
given in (35), for n — > 00. 

Tab. II shows the slope parameter d n (44) for 7 < n < 20, evaluated for 7 = 1 in various 
approximations; in the left part for R type orbits (odd n) and in the right part for L type orbits 
(even n). d^ 1 in columns 3 and 7 are the non-perturbative analytical results from (B42), d s ^ 
in column 2 represents the perturbative semi-analytical result (CIO) for the R orbits, and d™ um 
in columns 5 and 9 are the numerical results [10]. Columns 4 and 8 contain d™ um * obtained 
numerically from solving the equations of motion (3) and (4) for the periodic orbits with using the 
FA initial conditions (B12) and (B19) at the top turning point (Bl), and Eqs. (19) at t = T for the 
monodromy matrix elements. This approximation is in good agreement with the full numerical 
results for large enough n, the better the larger n, as seen from comparison of the 4th and 5th (and 




(43) 



d 



n 



trM ±fl , £ (e = 1) -2 



(44) 
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the last two) columns in Tab. II. The bifurcation energies for n > 12 were taken analytically from 
Tab. I. For smaller n, they were obtained by numerically solving equation trM^A(e n ) = 2 with a 
precision better than |trMj_yi(e n ) — 2| < 10~ 9 . As seen from this Table, one has good agreement 
of the asymptotic behavior of d n of the perturbative d s ^ and even better analytical results d™ 
as compared with these numerical calculations. It should be noted also that the slope parameter 
(44) of the perturbative approach (CIO) within the FA, see (CI), even without the correction (C4) 
to the periodic orbit UA(t), is in rather good agreement with the numerical results presented in 
Tab. II, especially for asymptotically large n, with a precision better than 5%. However, the second 
correction in (CIO) above the FA improves essentially the slope parameter (44) in this asymptotic 
region. As noted above, the asymptotic values of the perturbative d™ and the non-perturbative 
d^ 1 , as well as the numerical FA result for d™ um *, all converge sufficiently rapidly to the asymptotic 
analytical number crl(X) = 6.18199717 given by Eq. (35), in line with the analytical convergence 
found above from (B46). 

Fig. 7 shows good agreement between the analytical (B42), semi-analytical (CIO) and numerical 
solving the GHH equations (3), and (4) for classical periodic orbits and (19) for the monodromy 
matrix with FA initial conditions for L 12 and R i3 as examples. Both these curves agree very 
well with the asymptotic analytical slopes 0^(7) within a rather wide interval of 7 even for not 
too large n of the orbits mentioned above. This comparison is improved with increasing n, the 
better the larger n, which gives a numerical confirmation of the analytical convergence of the 
tr Mj_R 5 i(e, 7) (B42) to the asymptotic crl(i) (35) at the barrier e = 1 for any 7. For larger 7, 
one needs larger n in order to obtain convergence of all the compared curves. 

V. SUMMARY AND CONCLUSIONS 

In this paper we have investigated the bifurcation cascades of the linear A orbit in a class of 
generalized Henon-Heiles (GHH) potentials. We were able to derive analytical expressions for the 
stability traces trM^(e) of the A orbit and tiM± R ^ L (e) of the R and L orbits bifurcating from it 
as functions of the energy, which are asymptotically valid for energies close to the saddle at e = 1, 
i.e., in the limit where the bifurcations energies e n approach the saddle: e n — > 1. Our results 
confirm analytically the empirical numerical properties of the "Henon-Heiles fans" that are formed 
by the asymptotically linear intersection of the functions tr M± R ^ L (e) at e = 1, as given in Eq. (34). 
We found good agreement of our alternative non-perturbative and perturbative asymptotic results 
for tr M^^ i(e) with the numerical results. As a bonus, we have also obtained asymptotically exact 
expressions for the bifurcation energies e n of the A orbit in the GHH system, given in Eq. (33). 
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Our results can be interpreted in the sense that the non-integrable, chaotic GHH Hamiltonian 
becomes approximately integrable locally at the barrier, i.e., for e = 1. 

Both our approaches may be useful, also for more general Hamiltonians, for semiclassical cal- 
culations of the Gutzwiller trace formula for the level density [4], extended to bifurcation cascades 
with the help of suitable normal forms and corresponding uniform approximations [3, 29]. A nor- 
mal form with uniform approximation for two successive pitchfork bifurcations has been derived 
and successfully applied to the HH system in [12]. In future research, we hope to generalize the 
normal form theory to infinitely dense bifurcation sequences with the help of the results of [12, 13] 
and the theory of Fedoryuk [31, 32]. Hereby the "HH fan" phenomenon for the stability traces 
might be useful. 
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Appendix A: ASYMPTOTIC EVALUATION OF trM_ LA (e) FOR e -► 1 

To obtain the stability matrix M±a for the orbit A, we have to solve the linearized equation of 

motion (15) for small perturbations around the orbit in the perpendicular direction. Since the A 

orbit moves along the y axis, we have q = x, p = x and (15) becomes (3) which is already linear in 

x. We thus find from the non-periodic solutions of (3) with small initial values x = x(t = 0), 

xo = x(t = 0). Let us denote these solutions by x(t; xo, xq). The elements of M±a (we omit the 

subscript "_LA" for simplicity) are then given by 

x(Ta\Xo,0) , x(TA;0,Xn) 

M qq = lim 1 A ' °' ; , M qp = lim 1 A \ 0) , Al 

M pq = lim *( r **°'°) , Mpp = lim Mo) . (A2) 

Pq zo-0 x *o->0 x 

We could not find exact analytical solutions of (3) using the exact function yA(t) (6) for the 
A orbit, for which (3) becomes the Lame equation. Only at the bifurcation energies e n , one of 
its solutions is a periodic Lame function which has known expansions [22]. For the non-periodic 
solutions, no expansions could be found in the literature. We can, however, solve (3) if we instead of 
the exact yAif) use the approximation &U(£) given in (27), which becomes exact in the asymptotic 
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limit e — > 1, and for which (3) can be reduced to the Legendre equation as shown below. We 
proceed separately for the two time intervals < t < t 2 and t 2 < t < Ta, as specified after (25): 

a) < t < t 2 : Solve the equation 



x 1 (t) + [l + 2iY 1 (t)]x 1 (t) =0, 



with the initial conditions 



xi (0) = x = , ii(0) = x — > 0, 



and obtain xi^)- 

b) £2 < t < T\\ Solve the equation 



x 2 (t) + [l + 2 1 Y 2 (t)]x 2 (t) = 0, 



with the initial conditions 



(A3) 



(A4) 



(A5) 



£2(^2) = xi(i 2 ) , i 2 (t2) = xi(t 2 ) , 

and obtain x 2 {Ta)- 

To do so, we transform equations (A3), (A5) by defining the following variables: 

z\ — z , z 2 = z — 2 K(/c) . 

Then, the equations (A3), (A5) can be written compactly as: 

d 2 



dz 2 



Xi(zi) + [B + Ata,nh.\z i )}x i (z i ) =0, (i = 1, 2) 



where 



B = 



6(1 + 271/1) A _ 127(3/3-3/!) 



(2/3 - J/i) 

We next go over to the new variables 



(2/3 - y\) 



Si = tanh(zj) . (i = 1, 2) 



Then (A8) is transformed into the Legendre equation: 



(1 - s^) ^Xi(si) - 2si ^Xi(si) 4 



d 

d? 



1/(1/ + 1) - 



1 - s 



(A6) 



(A7) 



(A8) 



(A9) 



(A10) 



Xi{ Si ) = 0, (i = l,2) (All) 



with 



j u = iv / A+~B, 1/ = (-1 +iV4A - l)/2 



(A12) 
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The Legendre equation (All) has the solution 

xfa) = CuP^Si) + C 2i QZ( Si ) , (2 = 1,2), (A13) 

where Pjf(s) and Q£(s) are the associated Legendre functions of first and second kind, respectively, 
with real argument — 1 < s < +1 (see [19]). The initial conditions (A6) for 0:2(52) have the form 

x 2 (-s K )=x 1 (s K ), (—T—) = Hi^) ' ( A14 ) 

where 

s K = tanh K(k) = tanh(a K T A /2) . (A15) 

Solution (A13) of equation (All) for % = 1 with the initial conditions (A4) yields the following 
expressions for the coefficients Cn and C 2i : 

C u = x D»(s F ) Q»{s F ) , C 21 = -x D»{s F ) P?(s F ) , (A16) 

with 

s F = tanh F( V ,k), DZ(s f )= [a K (l - s 2 F ) W^s F )] _1 . (A17) 
Here W^(s F ) is the Wronskian 

W*{ 8 ) = W{Qt{s),P»{s)} = C%(8)±P?(8) - Pp(8)±C%(8) . (A18) 

Analogously, we solve equation (All) for i = 2 with the initial conditions (A14) and obtain for 
C\ 2 and C22 the following expressions: 

Cn = , D t (s F r Q ^'^ M , (A19) 

Wu \Sk) 

C^DSM ^'^f^ . (A20) 

The coefficients etj here are: 

«i = Q£(-a*)^M - PS{s K )^Qt{-8 K ) , 

a 2 = Q£(a*)^QJJ(-**) - Q^s^^M , 
a 3 = P^8k)^{-8k) - P${-8k)^PS{8k) , 

a 4 = P*(- 8k )±-QP(8k) ~ Q1Z(8k)^PP(-8k) ■ (A21) 
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Using X2 in (A13) at t — T\ and (A19), (A20) for the coefficients C12, C22, we obtain the following 
expression for M pp defined in (A2): 



M PP = [W»{s F )W»{s K )\ 



-1 



aiQZM^PZM+atPZM^PtisF) 



+ a 3 Q^s F )^Q^s F )+a 4 Pi:(s F )^Q^s F ) 
To calculate M qq defined in (Al), we solve equations (A3), (A5) with the initial conditions 



(A22) 



xi (0) = x -> , xi(0) = x = 0, 



(A23) 



and then take into account the condition (A6). Using the same steps as for M pp , we obtain M qq 
in the following form: 



M 



Using the following explicit expression for the Wronskian (A18), 

1 r(i + 1/ + a*) 



(A24) 



we now find for the sum of M qq and M pp 



{s 2 - 1) r(i + v - fj.) ' 



(A25) 



(A26) 



Note that the energy dependence comes through the quantities /i, v given in (A12) and sk in 
(A15) via the turning points yi(e) given in (9). As must be expected, the result (A26) does not 
depend on the initial point yo- 

We recall that the result (A26) has been obtained using the approximation (27) for the function 
VA{t), which is based on the asymptotic expression (25) for the Jacobi elliptic function sn(z, k), 
valid in the limit k — > 1. We can therefore simplify the above result by taking asymptotic limits, 
valid for e — > 1, of the quantities appearing in (A26). Since we have omitted the next-to-leading 
correction to (25), it is consistent to keep only the leading asymptotic terms. (An evaluation of 
all next-to-leading order corrections would lead outside the scope of this paper.) 

Using the asymptotic forms of the Legendre functions through hypergeometric series (cf. [19], 
Eqs. 8.704, 8.705, and 8.737), we obtain for the leading term in (A26) the intermediate result 



tr M ±A (e) w 2 Re (e^^F^) , 



(A27) 
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where the function Fa (7) is defined by 

Here the period Ta and the quantities a K in (8), and /x, z/ given in (A12) still depend on the energy 
e. Now, for e — * 1, all quantities in (A27) except T^e) have finite limits, easily found from the 
limiting turning points y\ — > —1/2, y 2 — > 1, 2/3 — > 1. In particular, we get the limits: 

a«->l/2, //->2Vl + 27, i/^^(-l+V 4 87-l)- (e -> 1) (A29) 

The limit of Fa (7) will be denoted by Fa (7) and its limiting phase by $a(7) 

Fa (7) - Fa (7) = \F A {i)\e™*™ . (e - 1) (A30) 

Its modulus can be given analytically as 

~ = v/cosh(47ryi + 2 7 ) + cosh(7rV48 7 - 1) 

' M >l V2 sinh[27rvT+27l ' 1 ' 



We discuss only positive values of 7 here; for 7 < 1/48, the function cosh (71-^487 — 1) becomes 



equal to cos — 487). The phase $a(7) is defined through Eqs. (A28) and (A30); it turns out 
to be negative for all 7 > 0. 

Using the above limits, we finally get from (A27) the asymptotic expression for trMj_A(e) given 
in Eqs. (29) and (30) of Sec. III. 



Appendix B: ASYMPTOTIC EVALUATION OF trM ±RjL (e) FOR e 1 

As mentioned in Sec. Ill A 1, the stability matrix Mi of a periodic orbit in a two-dimensional 
system is found by linearization of the equations of motion in the phase-space variables £± = (q,p) 
transverse to the orbit. For the R and L orbits, which have curved shapes that are only known 
numerically above their bifurcation energies, we have no way of determining the variables (q,p) 
analytically. We are therefore forced to evaluate the diagonal elements of the full monodromy 
matrix M, in order to find trM^ = trM— 2 for these orbits. For their calculation, we exploit some 
symmetry relations which are valid when the starting point at t — of a periodic orbit is chosen 
to be the upper turning point in the direction of the A orbit, i.e., along the y axis: 

y(t = 0) = y max , x(0) = . (Bl) 



We first present these relations for the A orbit and then for the R and L orbits. 
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1. Symmetry relations for elements of monodromy matrix M 

a. Diagonal elements for A orbit 

For the straight-line librating orbit A, we have r\\ = y, f\\ = y, and hence we may apply 
immediately (24). For the calculation of the elements M xx and M xx , we note that the differential 
equations for X xx (t) and X xx (t) contained in (19) decouple for the A orbit. Writing them at 
the time t = T, where T is the period of the A orbit, they can be combined into the following 
second-order differential equations for M XX (T) and M XX (T) as functions of the variable T: 

d 2 



dT 2 

d 2 



M xx (T) + V xx (T)M xx (T)=0, (B2) 



, r2 M xx (T) + V XX (T) M XX (T) = V xxy (T) M XX (T) y A (T) , (B3) 

where the subscripts of V denote its corresponding (successive) partial derivatives. With the 
special choice of the starting point (Bl), which for the A orbit becomes t/a(0) = 2/2, see (6), 
we have ?/a(0) = Va(T) = and the two equations for M XX (T) and M XX (T) become identical. 
For solving them uniquely, two boundary conditions are sufficient. Since both quantities become 
unity at bifurcations, may we choose two successive periods T = T n = T(e n ) and T = 2T n at the 
bifurcation energy e = e n to impose the boundary condition 

M xx (T n ) = M xx (T n ) = 1, M xx (2T n ) = M xx (2T n ) = 1 . (B4) 

This ensures the uniqueness of the solutions, so that we obtain the result 

M£> = M&> , (B5) 

which holds at arbitrary periods T and hence at arbitrary energies e. 



b. Diagonal elements for R and L orbits 

For the R orbits born at the successive bifurcation energies e n , we have ry = x, f\\ = x at the 
starting point (Bl), while y is the coordinate perpendicular to the orbit and one may apply (24). 
To obtain the elements M yy and of the R orbits at the starting point (Bl), we may use the 
"frozen approximation" (FA) for the y motion of these orbits (cf. [9, 10]) which is taken to be that 
of the A orbit, yn(t) ~ yA(t), so that the starting point is at y ma x = 1/2- This corresponds strictly 
to the lowest order of the perturbation expansion in the small parameter e = e — e n . Then, the 
velocity v x of their x motion close to e = e n is proportional to — e n as given in (B12) below. 
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For the functions M yy (T) and M yy (T), equations analogous to (B2) and (B3) hold, but with the 
subscripts x,x and y, y exchanged and jjA replaced by xr, T now being the period of an R orbit; 
boundary conditions analogous to (B4) apply. Hence we can conclude that in the limit e — > 1, 
where e = e — e n becomes small, the following approximate symmetry relation holds for the R 
orbits: 

M yf « M vf ■ ( e "> 1) ( B6 ) 
For the L orbits, the situation is slightly more difficult: their upper turning point does not 
lie on the y axis, nor do they reach or leave their turning point in the x direction. However, 
the x coordinate at the turning point is proportional to y/e — e n close to their bifurcation energy 
e n . Furthermore, the coordinate system (x,y) can be rotated such that the L orbits move in the 
rotated x direction at their upper turning points, and the diagonal elements of M are not changed 
under this rotation. Thus, the relation (B6) is, to leading order in e = e — e„, also found to hold 
for the L orbits: 

Mg) « Mg> . (e - 1) (B7) 



c. Relations of diagonal to non-diagonal elements 

Other symmetry relations can be obtained by taking the variational (partial) derivatives of the 
energy conservation equation at t — T: 

H[x(T),y(T),x(T),y(T)] = E, (B8) 

with respect to the initial variables, e.g., y(0) and y(0). Differentiating (B8) in y(0) and y(0) and 
applying the definition of the monodromy matrix elements (18), one has 

V X M X y + Vy Myy + X M ± y + ijMyy = Vy , 

V X M X y + Vy Myy + XM ± y + ijMyy = ij . ( B f ) ) 

where 

dV dV 
V x = — = x(l + 2 1V ), V y = — =y(l-y)+ 1X \ (BIO) 

according to the GHH Hamiltonian (1). All coefficients in front of the monodromy matrix elements 
are taken at the periodic orbit under consideration: x = x po (T) = x po (0), y = y po (T) = y po (0), etc. 
From (B9) at the starting point (Bl) for the R orbit, which in the FA is ^#(0) = 2/2, xr(0) = 0, 
one finds with y^(0) = 

M yy = 1 - V 4 M yy = M ^ (BH) 



23 

where 

"x = i(0)»^^> ^2 = 1/2(1-1/2)^^^^, (B12) 

see (BIO). The results in (B12), as well as all approximate relations given below, are valid in the 
FA in the limit e, e n — > 1 (with e > e n ) and are correct to leading order in y/\ — e n . From this 
one obtains the two approximate symmetry relations 



Mg>«Mg\ Mg>«^Mg). (B13) 



"2 



The first relation follows from the identical differential equations for the functions Mx y (T) and 
Miy(T) at the turning point (Bl) of the R orbits: 

M iy (T) + [1 + 2 iyR {T)\ M± y (T) = -2 7 x R (T) M yy (T) , 

M ±y (T) + [1 + 2 iyR {T)\ M ±y (T) = -2 7 x R (T) M yy (T) , (B14) 

according to (B6), and their zero initial values at e = e n . The second symmetry relation in (B13) 
can be proved directly through their definitions (18), 



M. 



mi 



Sx{0,e) _ /5x(0,e)/S( 



(B15) 



M yx 6y(0,e) \5y(0,e)/5e, 

where we write explicitly the energy dependence of the trajectory {x(t, e), y(t,e)} po owing to the 
initial conditions besides of the time dependence considered above. By employing the condition 
at the R top point, we find 

dy R (T, e n 



y R (T(e), e) = = y R (T(e n ),e n ) + (e - e n ) 



y R T'(e n ) + 



(B16) 



de 

We used here the y equation of motion (4) in order to obtain the derivative in the denominator 
of the r.h.s. in (B15). For the derivative in the numerator we may use the FA near the saddle 
energy, 5x(0,e n )/5e = ±(0)T'(0), because the main energy dependence is coming through the 
period T(e) in the argument of x po (T, e). Finally, from (Bll) and (B13) one arrives at two other 
useful approximate symmetry relations 



M^l + M^l + ^Mg). (B17) 



V, 
'2 



In an analogous way, from (B9), one directly derives the following symmetry relations for the L 
orbits accounting for their different initial conditions at the top (turning) point, |/l(0) = xl(0) = 0, 
y L (0) = J/2, x L (0) = x 2 (cf. [10]), 

Mg>«l + ^>«l-]|Mg>, (B18) 
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where 



V 1 = x 2 (l + 2 1 y 2 ), x 2 ^ ^/ 3(i e +2 e ^ — , v 2 = y 2 {l-y 2 )+ 1 xl^ ^j 1 -^, (B19) 

where the FA has been used. 

Other symmetry relations between monodromy matrix elements can be obtained in a similar 
way. In particular, one obtains the following structure of the stability matrix for both R and L 
orbits, 

•H</// M,A ^f M yy M yy Tl 



M ^=U M *U ±1 M J ' (B20) 

where the upper sign holds for R and the lower for L orbits. 

All approximate symmetry relations (B6), (Bll), (B13), (B17) and (B18) and the structure 
(B20) of the stability matrix have been checked by explicit numerical calculations, solving (19) 
in the FA with the given starting conditions. They become the more accurate the closer the 
bifurcation energies e n are to the saddle energy e = 1. 

In conclusion, we need not calculate those three quarters of the matrix X(t) in which the indices 
x or y appear. The coupled differential equations for the remaining elements of X(t) are 

X xx (t) + [l + 2 iypo (t)]X xx (t) = -2 1Xpo {t)X yx {t), 
X yx (t) + [1 - 2y po (t)} X yx (t) = -2 7 x po (t) X xx (t) , 

Xyy(t) + [l-2y p0 (t)}Xyy(t) = ~ 2^ X pQ (t ) X xy (t) , 
X X y(t) + [1 + 2 lyp0 (t)} X X y{t) = -2 7 X p0 (t) Xyy{t) , (B2l) 

to be solved with the initial conditions 

X XX (0) = Xyy(0) = 1 , X XX (0) = Xyy(0) = , 

Xxy(0) = Xyx(0) = , X^(0) = Xy X (0) = . (B22) 

In the equations (B21), the functions y po (t) and x po {t) describe the y and x motion of the periodic 
R and L orbits, respectively, born at the bifurcations. 

By using the relations (24) for r\\ = x, f\\ = x and (B6), (B7), one has the stability matrix trace 
of Mi for the R and L orbits, 

trM_L w 2 [X XX {T A ) + X yy {T A )\ -2=2 {M xx + M yy ) -2 = 2M yy , (B23) 

where T4 is the period of A orbit, T4 = T^(e), taken in the FA at the bifurcation energy, e = e n , 
= Xij(TA) (M xx = l). 



2. Analytical asymptotic expressions for tr Mxfl,L (e) 
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To solve the system of equations (B21), we have to specify the functions x po (t). In the asymp- 
totic limit e — > 1, we can use the FA in which y po (t) ~ 2m(0- The stability equation for the R and 
L orbits then is 

Xpo(t) + [1 + 2 iyA (t)} x po (t) = , (po = R, L) (B24) 

with the initial conditions xr(0) = 0, xr(0) = v x and xl(0) = x 2 , xl(0) = 0. As discussed in 
[10], (B24) with the exact yA(t) given in (6) is the Lame equation, whose periodic solutions are 
the periodic Lame functions. However, for e — > 1 we may replace the sn function in (6) by its 
asymptotic form given in (25): 

y A {t) nyi + (y 2 -yi)s 2 (t), s(t) = tanh [a K t - K(k)] , (B25) 

and transform the equation (B24) to the Legendre equation (All), replacing s, — > s and Xj(sj) — > 
Xp„(s), with s(t) given in (B25). We then obtain the x po (i) in terms of the Legendre functions as 

x R (t) = v x $+ (s K , s)/(a K W + ) , x L (t) = x 2 *+(s, -s^)/W+, (B26) 

where 

*±(s,s 1 )=Q%(s)Pg{s 1 )-Pg(s)Q%(s 1 ), (B27) 

g ^ ±(s) ^ ±(Sl) " P " ±±(s) o^ g " ±±(Sl) 
for the case of low plus index (minus will be used below). Here Q%±(z) and P£±(z) are the 
same Legendre's functions, as in Sect. A, sk is given by (A15), respectively. The constants W± 
independent of s is related to the Wronskian (A18), (A25) by 

T(l + v± + (j,±) 

r (i + v± - n±) 



*±(s,si) = 



[s\ - 1) (B28) 



w ± = ( S 2 - 1) wgw = ;:;:; ± "^ (b 29 ) 



with 



H ± = i*/A ± + B ± , i/± = i(-l+i v /4A^T), (B30) 

A+ = 127 (.m-yj B+ = 6(1 + 2-,,,) 

(j/3 - yi) (j/3 - j/i) 

, _ lafc-w) B _ = f-**) (B32) 



(ys-yi) ' (ys-yi) 

The solutions (B25) and (B26) are approximately periodic, the better the closer to the barrier 
energy. Note that in their derivations, we found more conveniently to use the initial conditions at 
t = Ta for R and t = for L orbits. All energy-dependent quantities, /i± and v± given in (B30), 
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k and a K in (8), as well as sk in (A15), are taken at the bifurcation energy e = e n like Ta in this 
approximation. 

For calculation of the trMj^z, (B23) through the symmetry relations (B17) and (B18), one 
has to derive the non-diagonal monodromy matrix elements M yx = X yx (TA_) and M xy = X xy (Tjy). 
Neglecting the right-hand sides of the first and third equations in (B21) and substituting their 
solutions 

Xg(a) = s K )/W +) X y ° y \s) = -s K )/W., (B33) 

into its second and forth equations, where the right-hand sides already contain small v x oc -Je — e n 
(B12) and x 2 oc -Je — e n of (B19) near the barrier, one finds for the solutions of the last two 
equations for X yx (t) and X xy (t), up to higher-order terms in the parameter a/1 — e n , 

X yx (t) = ^ r [ dt 1 x po (t 1 )^(s,s 1 )X^(s 1 ), (B34) 
a K W- Jo 

X xy (t) = -^=- [ t dt l x po (t 1 )<!> + (s,s l )X y y \s 1 ), (B35) 
a K W + Jo 

see (B28) for ^±(s, si) and (B27) for <E>±(s, si), with the same relation of t and t\ to s and si 
through s(t), see (B25), as explained above. In these derivations, we used the same transformation 
of equations of system (B21) to the Legendre form (All) via s = tanh(z) like above. 

With the help of (B34) for M yx = X yx (T A ), (B17) for M yy of R, and (B35) for M xy = X xy {T A ), 
(B18) for M yy of L orbits, and the periodic-orbit expressions (B26), by using the new variable s(t) 
of (B25) in (B34) and (B35), one obtains 

4 7 (e - e n ) f SK ds 



trM ±R (e) = 2 _ _ 2 / g $ + ( s *> s ) $ - ( s *> s ) ^+( s ' - s *)> 

3a3 i-sx 1 - s 

Ary (p _ p ) r s K Ag 

trM ±L (e) = 2 M _ J^ 2 / - tt + (s, -s K ) $ + (s K ,s) tf_(s,- SAr ). (B36) 

All factors, except for e — e n , on right of (B36) can be considered at the bifurcation energy e = e n 
for e n — > e — > 1. We neglected here, as in previous subsections of this Appendix, corrections of 
higher order in the small quantity y/l — e n . 

The integrals in (B36) can be taken analytically by simplifying their integrands with the approx- 
imation for the functions <3>_(sx,s) (B27) and ^-(s, — sk) (B28) with indices "— ", f-(e) — > —4 
and //-(e) — > —2 valid well in the limit e — > 1, see (B30) and [19], 

•-<*. *> - eo5?4) p "< s) - *-<«■ - s -> * "StS) p "W' < B37 > 



where the indices "— " appear only through a smooth energy-dependent coefficient, 

r(-//--i/--i)Q/- + //_ + !) 
r (n_ - 1/_ - 1) - + 1) ' 



(B38) 
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and T(e) tends to 120 smoothly at e — > 1, see (B38) and (B30). Notice that the contribution of 
the correction [19] to this approximation is negligibly small for the calculation of this integral, 



being of relative order y/1 — e n . Therefore, within the approximation (B37), the integrals over s 
in (B36) are reduced to the sum of several standard indefinite integrals of the products of two 
Legendre's functions with weight s and indices v = v + and fi = fi + of (B30) [33], 

dssCt{s)Zt{s) = K%(s) = Zj(s) -N 2 +r u (s)£» +1 (s)] 

- X 3 £: +1 (s)£» +l (s), (B39) 

where £^ and are any pair of the Legendre functions from the set P{f, Q%, 

Kl " M^T) ' * 2 " M^T) ' Ks " 2v{y + \) ■ (B40) 

The strong energy dependence of trMj_(e) (B23) is coming through the W_(e) or 1 — s 2 K {e) which 
tend both to zero for e — > 1 via the approximate key relations 



W4e) « - 35 ^gV r? , 1-^)"^, (B41) 

see (B29) for W_, (B30) for and z/_ [19]. The key point of our transformations is to remove 
indeterminacy zero by zero by identical cancellation of the singular factor W^(e n ) near the saddle 
from the denominators and that of the functions (B37) with indices "— " in the numerators of the 
integrands. Then, another constant singular factor 1 — s 2 K can be taken off the integrals. Thus, 
after such simple algebraic transformations with help of (B37) and (B39), from (B36) one obtains 

tr M ±R , L (e) = 2 ± 7 ( RL (e n ) (e - e n ) (l - s 2 K ) ^ p (s K )QZ(s K )±Q»(-s K ) + V^(s K ) 

(B42) 



where 



, Ssxhie) jl/s K for R\ 
, / n 1 ( r V s K 288 ^ , 

^ = (120) JT^f * W^)' (B44) 

vSi") = <fM - K§(- SK ), (B45) 

at s = sj<- with the indices "+", regular in the considered limit, fi = fi + and z/ = z/ + . Note that 
the derivatives of the Legendre functions on the right of (B42) are approximately proportional to 
1/(1 — s 2 K ), according to the recurrence relations for the Legendre functions with indices "+" of 
(B30) [19], and therefore, the product of the factor 1 — s 2 K by the expression in figure brackets is a 
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smooth function of the energy e n near the saddle as well as W+, see (B29) and (B30). Therefore, 
the strongest energy dependence near the saddle is coming only from the coefficient Crl(c) (B43). 
By making use of asymptotic expressions (B44) for 61(e) and (B12) for V^(e) through (B43) for 
Crl in the limit e n — > 1, up to higher order terms in small parameter a/1 — e n , from (B42) we 
arrive at the result 

trM w (e) = 2 ± 7 24(e ~ ^ \v^(s K )Q^s K )^(-s K ) + V^{s K )P^s K ) 

1 — e„ yj/^ as 

(B46) 



x^P?(-sk)-V$(8 K ) 



Q$M^P?(-8 K )+PJ!(8 K )±Qtt-8 K ) 



Using the asymptotic forms of the Legendre functions in figure brackets through hypergeometric 
series (cf. [19], Eqs. 8.704, 8.705, and 8.737), like for the derivation of tr M ±A (e) (A27), we may 



expand the function sk in (B46) in the small parameter 1 — s^(e„) oc a/1 — e n , see (B41), in the 
limit e n — > 1. Up to higher terms of relative order a/1 — e„, we then obtain the asymptotic result 
for trM^ L (e) given in (34), correctly describing the "GHH fans", with the slope function Crl{i) 
given in (35). 



Appendix C: PERTURB ATIVE CALCULATION OF trM ±i?L (e) NEAR e = 1 

I. "Frozen approximation" (FA) for the periodic orbits 

Within the FA, we set y po (t) m y A (t) (cf. Sec. Ill C). In order to find M yy of (B23) for trM ± , 
we solve the system of equations (B21) for X yy (t) and X xy (t), with the initial conditions (B22), 
iteratively by exploiting the smallness of their r.h. sides. Substituting expansions (41) and (42) into 
these equations at zero and first order in e, respectively, and then solving them for the monodromy 
matrix element M yy = X^T^), one obtains 

M ra = l + (e-e n )M« 1 (e n ). (CI) 

The first term is given by M yy (e n ) = 1 at order zero of the perturbation scheme, see (Bll). For 
the coefficient Vi y l y l {e n ) in the linear term of (CI) for the R type orbits, one finds 

KvM) = -7 2 6i(e„) Oe n ), (C2) 

where 61(e) is given by (B44), 

= f_ K dss<S> + (s K ,s){QZ(s K ) [Q^s)V P ^{s)-P^s)V^{s)] 

- pf{s K ) [Qtts)V%(8) - P»{s) Vg?( 8 )] } , (C3) 
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X^jf is the matrix (B45). For L orbits, one has similar derivations. All quantities on the r.h.s. of 
(C3) are taken at the bifurcation energy e = e n . In these derivations, the double integrals were 
reduced to single integrals through simple algebraic transformations with the help of (B37) and 
(B39), canceling the singular multiplier W_, see (B41), from the denominator with that in the 
numerator functions (B37) in the integrand near the saddle, like in Appendix B2. 



2. Corrections to the FA 

The results (C1)-(C3) can be improved much beyond the FA by taking into account the next- 
order terms in the expansion (41) of y po (t). We then find more exact solutions to (4). For instance, 
for the R orbits we get xn(t) = ur x^\t) and yii(t) = yA(t) + e Vii\t), which obey the initial 
conditions y^(0) = and y^(0) = 0, whereby 

VR(t) = T^fT f T%[4 0) (*i)r<M*i,«) (C4) 



3a*W- J- SK 1 - si 

with the relations s = s(t) and si = s(ti) of (B25). By making use of this solution and the 
corresponding more exact expansion (42) for X yy (t) of the problem (B21) and (B22), one finds 
correction to tr M±n(e). For these more exact calculations, we have to extend (CI) to the complete 
solution for M yy , collecting all leading corrections of first order in e, 

M yy = 1 + (e - e n ) [MW(e n ) + M«(e n )] , (C5) 

where 



; (e n ) = -^=- [ A dt 2 X < S(s 2 )^(s K ,s 2 )y^\t 2 ) (C6) 
dnVy ~ Jo 

with (C4) for y R ^(t 2 ). After a change of the integration variable from t 2 to s 2 = s(t 2 ) of (B25) 
in (C6) and using expression (C4) for y^{t 2 ), we may use the same approximations (B37), with 
the help of (B33) for X y ° y \s 2 ), to perform analytically the integral in (C4) in terms of elementary 
functions. Finally, after canceling identically the singular factor W _ and another singular factor 
1 — s 2 K from the remaining integral, like in Appendix B2, see (B41), one arrives at 

KUen) = -7 h(e n ) Oe n ), b 2 (e n ) = ±T W + h(e n ) : (C7) 



where 



C^n) = f' K ds^r£- {Pz\s)[T Q {s) - T Q {s K )\ 

- Qf(s)lF P (s)-F P (s K )]}, (C8) 
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5760 



30s - 118s 3 + 210s 5 - 90s 7 - 15 



In 



+ 192s 4 JFp( s ) In 



1 + s 
1 - s 



Fp{s) 



192 



1 + 8 
1 - S 



(6s 4 - 8s 6 + 3s 8 ) . 



(C9) 



Taking into account both energy corrections in (C5) with (C2) and (C7), we transform (C5) 
for tr M±r into the asymptotic result 



tr M ±R (e) = 2 - 2(e - e n ) 7 2 h{e n ) 1^^) + lb 2 (e n )I^ 2 (e n ) 



r (i) 



(CIO) 



L(l-e n )V2_ 

Similar expression for the stability trace trMj_^(e) can easy obtained for the L orbits. As seen 
from (C7) and (B44), 6 2 (e n ) oc &i(e„) oc 1/(1 — e n ), and the other factors Z^i(e n ) (C2) and 
(e„) (C8) are smooth functions of e„, as confirmed by numerical integrations in (C2) and 
(C8). Therefore, both corrections in (CIO) are mainly proportional to (e — e n )/(l — e n ), i.e., linear 
in the (e — e n ). They are both finite in the barrier limit e — + 1 but numerically, the essential 
contribution to (CIO) is coming from the first FA correction while the second one (above FA) 
is much smaller. Note that the leading higher-order terms in the parameter e (40), originating 
from the next iterations in the perturbation scheme (41) and (42), can be estimated, in fact, as 



of higher order in y/1 — e n . Thus, the complete sum of energy-dependent corrections in (C10) 
has the same leading energy dependence oc (e — e n )/(l — e n ), up to higher-order terms in small 
parameter \/\ — e n , as in (34). The leading energy dependence of trMj_(e) is thus precisely that 
found explicitly in the non-perturbative result (B46) in Appendix B 2. 
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Figure 1: Trace of stability matrix Mi of orbit A and the orbits born at successive pitchfork bifurcations 
in the standard HH system (7 = 1), plotted versus the scaled energy e. From bottom to top: successively 
zoomed energy scale near e = 1 (from [9]). 
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Figure 2: Trace of the stability matrix of the orbits A (heavy line), B, C, and the orbits R2 m -i, 

(m > 3) born at successive pitchfork bifurcations of orbit A in the standard HH potential, plotted versus 

their individual periods T. AT is the asymptotic period of the curve tr M±a (Ta) for large Ta (from [9]). 



34 




Figure 3: The "Henon-Heiles fans". Trace of stability matrix of primitive A orbit (solid line) and the first 
four pairs of R orbits (dashed) and L orbits (dash-dotted) in the standard HH system, plotted versus 
scaled energy e; the latter forming two fans for the R and L orbits. 




Figure 4: Stability discriminant trM^A of the A orbit in the HH potential (7 = 1), plotted versus 
period Ta. Solid line: numerical result (as in Fig. 2, from [9]). Dotted line: analytical asymptotic result 
trMfJ(T A ,l) given in (31). 
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Figure 5: Upper panel: slope parameter crl of the "HH fans" plotted versus the potential parameter 7. 
Crosses: numerical values; solid line: the function crl{i) given in (35). Lower panel: excerpt for small 
values of 7. The dotted line gives the linear approximation to crl(i), with the slope given in (36), as 

fmmrl in a comir*! accural norfn lrKafnro annvna r*V» fill 




Figure 6: Stability traces trM_|_^,L — 2 as functions of the energy e at 7 = 1. Solid lines show the 
analytical expression (B42) for R n and L n orbits with n = 7 — 14. Dashed lines are the perturbative 
results (CIO) for a few R orbits as examples. 




Figure 7: Stability traces tr Mj_^l — 2 as functions of 7 for R and L orbits, respectively, evaluated at the 
barrier energy e = 1. Solid lines show the analytical expression (B42) for the orbits R13 and L12; dashed 
lines the asymptotic results (35) for Tcrl(t); dots are the perturbative results to (CIO); and crosses show 
the numerical results =F<i™ um * with the FA initial conditions as in Tab. II. The bifurcation energies e n {j) 
are obtained analytically through Eqs. (33). 
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n 


e* 


^n 


5 





96945 


52246 


81049 





96930 


90904 


6 





98682 


99363 


40510 





98670 


92353 


7 





99918 


81219 


03970 





99918 


78410 


8 





99964 


99405 


84051 





99964 


98 


9 





99997 


84203 


34217 





99997 


8390 


10 





99999 


06954 


44011 





99999 


06955 


11 





99999 


94264 


13919 





99999 


9424 


12 





99999 


97526 


85521 





99999 


97525 


13 





99999 


99847 


54120 





99999 


99847 5 


14 





99999 


99934 


26398 





99999 


99934 3 


15 





99999 


99995 


94766 





99999 


99996 046 


16 





99999 


99998 


25274 





99999 


99998 249 



Table I: Bifurcation energies in the standard HH potential (7=1). e* : asymptotic values, calculated 
from the analytical expressions (33) up to 15 digits with MATHEMATICA. e n : numerical values 
taken from [10]. 



n 


a n 


Jan 
u n 


jnum* 
a n 


r inum 
u n 


n 


jan 
u n 


jnurrt* 
u n 


r inum 
u n 


7 


4.7476 


5.5863 


6.1688 


6.1801 


8 


5.7796 


6.2661 


6.1803 


9 


5.9234 


6.0901 


6.1800 


6.1819 


10 


6.1209 


6.1951 


6.1820 


11 


6.1391 


6.1685 


6.1817 


6.1820 


12 


6.1731 


6.1841 


6.1897 


13 


6.1750 


6.1801 


6.1819 


6.1837 


14 


6.1807 


6.1823 




15 


6.1808 


6.1817 


6.1820 




16 


6.1818 


6.1821 




17 


6.1818 


6.1820 


6.1820 




18 


6.1820 


6.1820 




19 


6.1820 


6.1820 


6.1820 




20 


6.1820 


6.1820 





Table II: The slope parameters d s ^ of the semi-analytical (C10) and c^ n of the analytical (B42) 
expressions vs the numerical values d"" m * for solving (3), (4), (19) within the FA for the initial 
conditions at the top point (Bl) of the periodic R n (left) and L n (right) orbits, and d™ m is exactly 
full numerical results [10] (7 = 1 in all cases). 



